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ABSTRACT 

The class of tidal features around galaxies known variously as "shells" or "umbrel- 
las" comprises debris that has arisen from high-mass-ratio mergers with low impact 
parameter; the nearly radial orbits of the debris give rise to a unique morphology, a 
universal density profile, and a tight correlation between positions and velocities of 
the material. As such they are accessible to analytical treatment, and can provide 
a relatively clean system for probing the gravitational potential of the host galaxy. 
In this work we present a simple analytical model that describes the density profile, 
phase-space distribution, and geometry of a shell, and whose parameters are directly 
related to physical characteristics of the interacting galaxies. The model makes three 
assumptions: that their orbit is radial, that the potential of the host is spherical at 
the shell radii, and that the satellite galaxy had a Maxwellian velocity distribution. 
We quantify the error introduced by the first two assumptions and show that selecting 
shells by their appearance on the sky is a sufficient basis to assume that these simpli- 
fications are valid. We further demonstrate that (1) given only an image of a shell, the 
radial gravitational force at the shell edge and the phase-space density of the satellite 
are jointly constrained, (2) that combining the image with measurements of cither 
point line-of-sight velocities or integrated spectra will yield an independent estimate 
of the gravitational force at a shell, and (3) that an independent measurement of this 
force is obtained for each shell observed around a given galaxy, potentially enabling a 
determination of the galactic mass distribution. 
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1 INTRODUCTION 

Recently, large-scale sky surveys and deep follow-up images 
have been used to discover a wealth of tidal debris around 
our Galaxy and others nearby (e.g., Ibata et al. 2001; Mc- 
Connachie et al. 2009; Trujillo et al. 2009; Martmez-Delgado 
ct al. 2010; Radburn-Smith et al. 2011). This tidal de- 
bris comes in many shapes and sizes, from the huge tidal 
arms created by interacting, roughly equal-mass galaxies to 
the fainter features observed around nearby galaxies. These 
smaller-scale features are thought to come from interactions 
of the large galaxy with much smaller satellites, which are 
called minor mergers. Evidence that minor mergers occur 
in nature is an important link to our cosmological history, 
since cosmological simulations of dark matter indicate that 
about half the mass in the Milky Way's outer regions was 
accreted in this way (for example, Maciejewski et al. 2011; 
Wang et al. 2011). Minor mergers are also useful as a way 
to constrain the shape and mass of the large galaxy, since 
tidally stripped material from the smaller satellite galaxy 
behaves as test particles in the relatively undisturbed po- 
tential of the larger host galaxy (Ibata et al. 2001; Helmi 
2004; Johnston et al. 2005; Eyre & Binney 2009; Law & 



Majewski 2010). The remnants of these mergers give us a 
way to measure the characteristics of the dark components 
of galaxies, which are predicted with great accuracy by cos- 
mological models and simulations. 

Some minor mergers create patterns of tidal debris that 
look like shells or umbrellas. The first such debris was iden- 
tified around elliptical galaxies by Malin & Carter (1983); 
these galaxies were called "shell galaxies" because of these 
distinctive features. Hernquist & Quinn (1987, 1988, 1989) 
showed that the shells were probably created by a minor 
merger on a nearly radial orbit. This explains the alternate 
spacing of the shells on either side of the host galaxy since 
they are formed as material from the satellite piles up ap- 
proximately at turning points: the satellite initially had a 
distribution of energies that is reflected in the different radii 
of the shells. More recently, similar features have been dis- 
covered around nearby disk galaxies (Ibata et al. 2001; Mc- 
Connachie et al. 2009; Martmez-Delgado et al. 2010). The 
vast improvements in imaging since shell galaxies were first 
identified, and the relative proximity of these objects, have 
revealed more of their structure than had previously been 
observed. In some cases (Fardal et al. 2012; Romanowsky 
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et al. 2012), the objects are even close or bright enough that 
velocity information could be obtained. 

Shells from nearly radial mergers are particularly spe- 
cial because there is a direct correlation between the kine- 
matic properties of debris and its location relative to the 
host galaxy. Thanks to the near-symmetry of the encounter, 
the system can be considered in a two-dimensional projec- 
tion (r,v r ) of the full six-dimensional phase space without 
much loss of information. In this two-dimensional projection, 
the initially cold satellite material, once unbound, forms a 
thin stream that winds through phase space, so that for any 
spatial location r there are a small, finite number of streams 
with different characteristic v r - Merrifield & Kuijken (1998) 
pointed out that this correlation could be used to estimate 
the mass of the galaxy hosting a shell, if the line-of-sight ve- 
locity of the material could be measured at different points 
in the shell. 

The dynamics governing the formation and shape of 
this stream are closely related to earlier work on spherically 
symmetric secondary infall of matter accreting onto dark 
matter halos. As shown in Fillmore & Goldreich (1984) and 
Bertschinger (1985) radial accretion of gravitating, cold, col- 
lisionless matter forms a series of infinite-density peaks at 
successive radii, known as caustics. Mohayaee & Shandarin 
(2006) further showed that for warm matter with a finite ve- 
locity dispersion, the peaks take on a finite width and height, 
but although they are no longer caustics in the mathemati- 
cally rigorous sense they retain many of the same properties. 
This work was motivated by cosmological simulations and 
so considered the radial collapse of spherically distributed 
matter, but Helmi & White (1999) showed that caustics 
are also produced when a small, initially self-bound satel- 
lite falls into the (assumedly) static potential of a larger 
host galaxy. In fact, caustics are a universal product of dy- 
namical systems with turning points, regardless of the type 
of symmetry (Tremaine 1999; Hogan 2001). This work has 
been recently confirmed with state-of-the-art computer sim- 
ulations using realistic galactic potentials by Vogelsberger 
et al. (2008), and systems resembling shell galaxies are pro- 
duced in cosmological simulations with semi-analytic stellar 
components (Cooper et al. 2011). Sanderson &: Bertschinger 
(2010) demonstrated the connection between the shape of 
the density profile of the caustic and the physical character- 
istics of the interacting galaxies. For these reasons, in this 
work we will refer to shells as "tidal caustics" , a term which 
emphasizes their high degree of symmetry and the correla- 
tion of positions and velocities. 

Recently, Ebrova et al. (2012) have discussed the possi- 
bility of using integrated-light spectra to measure the grav- 
itational potential in shell galaxies by obtaining the line-of- 
sight velocity profiles of the shell debris. Their results are 
similar to some of those given in this work, with a few im- 
portant differences. First, they consider only spherical shells 
from satellites on perfectly radial orbits, not the effect of 
relaxing these assumptions to include angular momentum, 
potential flattening, or projection effects. In this work we 
explore how all three of these things affect the line-of sight 
velocity distribution. Furthermore, the line-of-sight veloc- 
ity profiles derived by Ebrova et al. do not account for the 
nonzero thickness of the shell. This causes the peaks of the 
profile in their model to be thinner than and slightly off- 
set from the peaks in their simulated shells; we show in this 



work that including the shell thickness solves both problems. 
Finally, they do not discuss how the surface-brightness pro- 
file of the shell is related to the kinematic profile, whereas 
we show that the two can be used together to simplify the 
fitting process by constraining the shell's geometry indepen- 
dently prior to modelling its kinematics. 

In this paper we present a simple analytical model for 
the density and phase space distributions of tidal caustics 
and relate the parameters of the model to characteristics 
of the interacting galaxies (Sections 3 and 4), specifically 
the radial component of the gravitational force g s produced 
by the host galaxy at the radius of a tidal caustic and the 
initial phase space density of the satellite galaxy. We dis- 
cuss how the geometry of the caustic can be modeled in or- 
der to calculate projections of the two analytic distributions 
(Section 5) and then present three examples that combine 
the phase-space and geometric models to transform into ob- 
servable space: surface-brightness profiles (Section 5), point 
measurements of the line-of-sight velocity (Section 6), and 
velocity profiles for spectra (Section 7). In Section 5 we com- 
pare the model qualitatively to simulations; in Section 6 we 
also demonstrate that g s can be recovered within 50 percent 
or better for most cases. In Section 7 we discuss the advan- 
tages of integrated-light spectra over point sources and in 
Section 8 we summarize. 



2 SIMULATIONS OF TIDAL CAUSTICS 

To check that the analytical model represents the salient 
features of caustics, we will compare it with a set of simula- 
tions intended to gradually relax the assumptions of spher- 
ical symmetry and radial orbits that are fundamental to 
its derivation. The simulations are labeled A through D in 
order of increasing realism, for comparison with the ana- 
lytical model; simulations B and C are actually two series 
that gradually increase the angular momentum in two dif- 
ferent potentials. Table 1 summarizes the parameters of all 
the simulations. 

In all four types of simulations we take a coordinate 
system in which z is the direction along the line of sight, 
centered on the host galaxy. In this case the x-y plane is ef- 
fectively the plane of the sky, and will be used to represent 
the view of the tidal debris that we would observe. For this 
work we assume that the angular size of the galaxy is small 
enough that we can ignore the slight difference between the 
z coordinate and the radial line of sight (spherical projection 
effects). Then, for simulations A through C, we choose the 
initial position and velocity vectors for the satellite galaxy 
to lie mainly in the sky plane, so that the resulting shells 
are symmetric with respect to z = 0. This is the orienta- 
tion, with respect to our line of sight, in which the shell 
edges appear sharpest and the kinematic signature is least 
complicated. In simulation D we relax this assumption. 

Simulation A (Figure 1, top left) is the base case that 
fulfills all the assumptions in the derivation: a simulation 
of the disruption of a self-bound satellite in a static, spher- 
ical potential on a purely radial orbit (Lcom = 0). The 
series of simulations labeled B uses the same potential but 
gradually increases the orbital angular momentum, param- 
eterized by the angle 9 between the initial position and 
velocity vectors of the center of mass. = corresponds 
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to a radial orbit; 8 = 90° is a circular orbit. Using this 
parameter, the angular momentum is Loom = -Lcirc sin #, 
with Lcirc = rcoMfc^coM). An example with a moderate 
amount of angular momentum, 6 = 20°, is shown in the top 
right panel of Figure 1. 

Series C uses a highly flattened, axisymmetric poten- 
tial and sets the initial center-of-mass velocity vector of the 
satellite at 45 degrees to the plane of symmetry (the x-y 
plane) for maximum sensitivity to the nonsphericity of the 
potential. The amount of orbital angular momentum is ad- 
justed as in series B, again parameterized by the angle 9. The 
result with 6 — is shown in the lower left panel of Figure 
1. This set of simulations is intended to exaggerate the ef- 
fect of departures from spherical symmetry, thus stretching 
the capabilities of the model to the limit. Indeed, we find in 
this case that caustics selected in (r,v r ), coordinates which 
no longer reflect the symmetry of the system, include both 
material passing through a turning point in R and material 
reaching a turning point in z. This further complicates the 
kinematic signatures for this case, but the projected distri- 
bution in the sky plane still features shells. 

Simulation D (Figure 1, bottom right) differs from the 
other simulations in that it was not a product of systematic 
experiments, but is instead an N-body model of a real set 
of caustics (the Andromeda Giant Stream and associated 
shells) in a realistic, though still static, galactic potential (a 
mass model of M31). The best-fit model obtained by Fardal 
et al. (2007) includes two caustics produced by the nearly- 
radial encounter of an intermediate-mass satellite galaxy. 
This simulation is included in lieu of actual data so that 
we can compare the full phase-space information inferred 
by the N-body model with our analytic expressions, while 
emphasizing that these more realistic conditions can still 
give rise to suitably symmetric caustics. 



3 A UNIVERSAL DENSITY PROFILE 

In previous work (Sanderson & Bertschinger 2010) we de- 
rived a simple analytical form for a one-dimensional caustic 
formed from a system with initial random velocities drawn 
from a Maxwellian distribution (Figure 2, left panel). In the 
case where this system is a small, gravitationally self-bound 
satellite galaxy falling radially into the center of a static host 
galaxy, the one-dimensional density profile of the caustic, as 
a function of the galactocentric radius r, can be written in 
terms of four physical parameters S r , r s , k, and /q: 



p(r) 



fo 



,r/isi 



B 



(r-r s ) 2 



4<5 2 



(1) 



where B is a piecewise combination of modified Bessel func- 
tions of the first kind: 



B{u) = 



1-1/ a{u) +l l / i {u) 

T_i/ 4 (m) -Zi/ 4 (it) 



r < r s 
r > r s 



(2) 



In Equation (1), 5 r is the characteristic width of the 
caustic surface, which depends on the initial velocity disper- 
sion of the dwarf galaxy that created the caustic. A perfectly 
"cold" distribution, in which all particles had the same ini- 
tial energy, would give rise to a true caustic with zero width 
and infinite density. We refer to caustics in a less rigorous 
sense: they have finite density and width but still exhibit a 



large local density enhancement thanks to the same dynam- 
ical processes. The width S r is related, by conservation of 
phase space volume, to the phase space distribution of the 
satellite, but the relationship is not straightforward. The vol- 
ume of phase space that is conserved is that of the material 
after it has been stripped from the satellite at pericenter. 
For such a radial orbit this volume can be roughly related 
to the pre-stripping phase space distribution of the satellite 
material via impulse approximation methods, but this in- 
volves further assumptions about the nature of the satellite. 
Likewise, the methods described in Helmi & White (1999) 
could in principle be used to relate the phase space volume 
to the width 5 r , but in practice many other factors (includ- 
ing projection effects, angular momentum, and departures 
from spherical symmetry in the potential) can increase 8 r , 
further complicating the relationship with the nature of the 
original satellite and increasing the number of assumptions 
necessary to obtain what amounts to a lower limit on the 
phase space density. For this work, we wish to remain as 
model-independent as possible with respect to the interact- 
ing galaxies, so we will treat S r as a free parameter in the 
remainder of this paper. 

The radius of the caustic surface is denoted r s . It is close 
to, but not equal to, the radius of peak density, r max , be- 
cause of the thickness of the caustic, induced by the nonzero 
velocity dispersion of the material. The peak location r max 
and the peak density p max may be determined by solving 
dpi dr = numerically (necessary because of the Bessel func- 
tions). The peak radius r max is related to the caustic radius 
r s by 



r s - 0.765 5 r , 



(3) 



and the density at this location is 



1.021 /ot/-l. (4) 
The curvature k measures the shape of the stream near 



r s in the (r, v r ) projection: 



dv% 



(5) 



Because by definition the caustic surface is located where the 
distribution is vertical in the projection shown in the right 
panel of Figure 2, near the caustic its shape may always be 
approximated by a quadratic function: 



(6) 



where v s is the radial velocity at the caustic surface. Al- 
though all the material forming the caustic is near apoc- 
enter, the net radial velocity at the caustic surface is not 
necessarily zero because the spatial overlap of the orbits of 
individual stars in the stream depends on the distribution 
of their orbital periods, which in turn depends on the dis- 
tribution of their energies. 

We can relate k to the gravitational force at r s by differ- 
entiating the energy equation evaluated at the caustic sur- 
face. For a purely radial orbit in a spherical galactic potential 
$(r), we have 



dr 



dv r 
dr 



dr 



9(r) (7) 



for the first derivative. If we invert the derivative in this 
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Table 1. Information about the four simulations whose results arc presented in this work. See text for additional discussion. Expressions 
for the potential types in the table are given in Appendix A. 



Name 



Potential type 



Potential parameters 



Satellite mass Initial tqqm (kpc) 
(xlO 8 M Q ) 



Initial vrjQM (km s 1 ) 
v c = » c (r COM ) 



A 

Series B 
Series C 
D 



Spherical Isochrone 
Spherical Isochrone 

"Flattened" Isochrone 

see Geehan et al. (2006) 



M = 2.7 x 10 12 Mq, b = 8.0 kpc 2.2 
M = 2.7 x 10 12 Mq, b = 8.0 kpc 2.2 



M = 2.7 x 10 12 M e ,b = 8.0 kpc, 
q = 0.7 

M M3 i(< 125 kpc) = 7.3 x 
lO 1:L Af0, rjjaio = 7.63 kpc; see 
Table 2 of Fardal et al. (2007) for 
complete parameters 



2.2 



22 
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(-40,0,0) (v c cos 6, v c sind, 0), 

9 6 {10° ...90°} 
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(-34.75, 19.37, -13.99) (67.34, -26.12, 13.50) 
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Figure 1. Views of the sky (x-y) plane for the four types of simulations used in this work. In each simulation, caustics selected in r — v r 
space (shown in different colors) reflect the shell features seen in the plane of the sky. Simulations A-C have the same mass ratio of about 
10 — 4 ; simulation D has a mass ratio of about 3 X 10~ 3 . 
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r s -105r r s -5i5r r s -2Sr r s r s +2Sr r s +5Sr 



r 

Figure 2. Left: Universal caustic form for radial infall from a population with initial Gaussian velocity dispersion. r 3 and 8 r are defined 
and formulae for the peak radius and density are given in the text. Right: Phase space distribution taken from simulation D. Near the 
peak of one of the caustics (black), the phase space distribution follows the curve given by Equation (6) (dashed line), where v s is the 
radial velocity of particles at the caustic surface, k is defined in the text. 



expression and differentiate with respect to v r , then express 
the first derivative in terms of g, we find that 



2<?(r s 



dg 



9(7 



dr 



(8) 



We can suppress the dependence on the tide, dg/dr, by us- 
ing Equation (8) to obtain a lower limit on k instead of 
an equality, under certain assumptions about the mass dis- 
tribution. It can be shown from the Poisson equation that 
dg/dr is negative for any spherical density distribution that 
satisfies the relation 



M{< r) 



p(r) < -p(r) 



(9) 



at the radius of interest. For a power-law density distribu- 
tion, par 1 , this condition is satisfied for all 7 < —1. Ob- 
servational evidence and simulations of dark halos both in- 
dicate that the density almost certainly falls off faster than 
this at most radii, especially in the outer regions where shells 
are usually found. For example, the log-slope of the density 
of the Milky- Way-like halos of the Aquarius project satis- 
fies this condition everywhere except in the innermost bin 
of the highest-resolution simulation (Springel et al. 2008) . If 
we assume dg/dr < 0, then the second term in Equation (8) 
is always positive, so that the quantity in brackets is always 
larger than 1. The Poisson equation also requires that g < 
at all r for a spherical potential, since there is no negative 
mass, implying that k is positive. Thus neglecting the term 
involving the derivative of g turns Equation (8) into a lower 
limit on n: 



k ^ 



1 

2|ff(r.)| 



(10) 



For most reasonable models of the halo the neglected term 
is much less than 1, given that v s is close to zero and the 
tidal forces are weak at large radius (where shells are most 
easily observed) . Neglecting the gravity gradient term gives 



(11) 



2|3(r s )r 



relating the curvature of the phase-space stream directly to 
the gravitational force exerted by the host galaxy. 



3.1 The effect of nonzero angular momentum and 
deviations from spherical symmetry 

If the satellite galaxy had some initial orbital angular mo- 
mentum L, the energy equation (7) is no longer correct, so 
measuring g by using the approximate relation (8) will have 
additional error. In this case, following the same derivation 
for k while including nonzero L (but still assuming a spher- 
ical potential) leads to the relation 



k — — - 



1 + 



rigs 



rlg a 



1 - 



2 + 



flf dr 
r s g s 



rigs 



(12) 



where g s = g(r s ). The overall l/2g dependence is preserved, 
and the first two terms in the square brackets correspond 
to those in Equation (8). However, there are now new terms 
that all depend on the dimcnsionlcss quantity 



rigs 



L 

^2^,2 ' 



(13) 



which compares the orbital angular momentum to that of 
a circular orbit at r 3 , since v c (r s ) = ^Jr s \g s \ is the circular 
velocity at r s . The minus sign is included in the definition 
to ensure a positive el, since g s is negative. When the orbit 
is purely radial, e L = 0; when it is circular, e L — 1. 

Solving Equation (12) for g gives a complicated depen- 
dence on the angular momentum that can be simplified by 
expanding in terms of el and its counterpart for the tidal 
force, 



£ t = 7, -r 



gi dr\ 



dg\ 

dr 



(14) 



and requiring that both be much less than 1. To first order 
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in both quantities, 



1 + e T - e L 1 + 



6kv, 



+ 0(e L e T ). (15) 



The leading-order change in re due to the angular momen- 
tum is at the L 2 level, so we expect that Equation (11) will 
remain a decent approximation even for satellites on fairly 
non-radial orbits. Furthermore, the corrections for tides and 
angular momentum have opposite signs, so at intermediate 
L we may expect help from competing error terms. 

In practice we will not be able to make even first-order 
corrections for either tides or angular momentum, since one 
requires independent measurement of the gravity gradient 
and the other requires knowing the initial angular momen- 
tum of the satellite. However, with numerical experiments 
we can test how sensitive the phase-space shape is to nonra- 
dial orbits. To do this, we obtain re from fits to the projected 
(r, v r ) distribution for shells in the simulations of Series B. 
We first estimate g using Equation (11), then progressively 
correct this estimate at first order according to Equation 
(15), taking first tides, then angular momentum into ac- 
count. The simulations in the series form four caustics at 
different radii, and we obtain an independent estimate from 
each. 

As seen in the left panel of Figure 3, the spherical ap- 
proximation still recovers the gravitational force to within 
about 20 percent, even for relatively large amounts of angu- 
lar momentum (up to half the maximum value). This is due 
to the quadratic dependence of e_L on L and the partial can- 
cellation of ex and . In some cases the combined first-order 
corrections for tides and angular momentum are sufficient to 
eliminate nearly all the error, while in others the correction 
is less good, especially at larger L. This is because the av- 
erage L of material in each shell is not necessarily equal to 
Lcom, which was used to calculate the correction. Since sl 
is quadratic in L this difference is magnified at larger L. 

We also see that the difference between a system that 
has less than 20 percent error using the radial approxima- 
tion (Figure 3, center panel) and one that has much more 
(Figure 3, right panel) coincides with the difference between 
a shell-like and a stream-like morphology. The example with 
more angular momentum still has a sharp edge at the apoc- 
enter of each caustic, but the morphology of the debris is 
otherwise more stream-like than shell-like. Deviations from 
spherical symmetry in the potential will populate a larger 
proportion of orbital planes per unit L thanks to precession, 
and this should further accentuate the differences between 
low-L systems and those with higher L. Projection effects 
will also increase these differences; the two systems shown 
here are at an ideal viewing angle where the shells look the 
sharpest. We conclude that identifying shells (consisting of a 
fan shape and a sharp outer edge) by eye as low-L systems is 
probably a sufficiently strict selection criterion for applying 
the radial-orbit approximation to obtain estimates of g s . 

Deviations from spherical symmetry in the potential 
can also affect the ability to recover g s . Out of the plane 
of symmetry, orbits in an axisymmetric potential will pro- 
cess, thickening the phase space stream perpendicular to the 
(r, v r ) plane. Futhermore in an axisymmetric potential, the 
spherical radius r and its conjugate velocity v r are not the 
coordinates in which the problem is separable, so there is no 
guarantee that the analysis developed here will still be ap- 
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Figure 4. As in the left panel of Figure 3, but for simulations in 
series C (a highly flattened potential). 



plicable to non-spherically-symmetric systems. To see how 
well we can do for a significantly flattened system, we repeat 
the process used to obtain Figure 3 on Series C. The results 
are shown in Figure 4. Equation (15) does not include the 
effect of precession, so the error prediction from the spheri- 
cal analysis does not fully correct the estimate, but we can 
still recover g within 20 percent for the cases shown. We 
conclude that although we expect galactic potentials to be 
flattened, the inaccuracy from assuming spherical symmetry 
is not large enough to prevent us from measuring g. 



4 A SIMPLE PHASE-SPACE DISTRIBUTION 

Caustics have a simple structure in phase space as well as a 
simple density profile, thanks to the fact that the material 
in them was initially compact in phase space. The preser- 
vation of this initial small phase-space volume means that 
stripped material will spread out in a thin stream through 
phase space; thus the phase space distribution at late times 
is nearly one-dimensional. In the case of the caustics formed 
by a nearly radial encounter, symmetry tells us that the 
stream flows mainly along two of the six available coordi- 
nate directions: galactocentric radius r and radial velocity 
v r . This is why, following Equation (6), we can construct 
a preliminary model of the caustic phase space distribution 
as a one-dimensional function of the six-dimensional phase 
space coordinates (x, v) by ignoring everything but these 
two coordinates: 



f(x, v) oc 5 [r s — r — hi(v r — v s ) 2 ] 



(16) 



Normalizing this distribution so that the integral over a shell 
is unity results in the expression 



f(r, Vr 



16r 5 s /2 n, 



5 [r s — r — re (v r — v s ) 2 ] 



(17) 



where f2 s is the solid angle spanned by the shell. The deriva- 
tion of this expression is given in Appendix B. 

Unlike the model of the density profile, this model for 
the phase space distribution does not self-consistently take 
the energy spread into account, since it will also result in 
a spread of the phase space distribution around this line. If 
Equation (17) is integrated over v r , one obtains the limit of 
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Figure 3. Left: the fractional error in recovering g remains small even at fairly large values of the orbital angular momentum (proportional 
to sin 8). The blue open circles show the relative error from using Equation (11) to estimate g. The red open circles include the first-order 
correction for tidal forces from Equation (15); the solid black points correct for both tides and angular momentum at first order using 
the same equation. Center and right: the tidal features at 8 = 50° (right), when the error on g B exceeds about 20 percent, still exhibit a 
sharp outer edge but are morphologically distinct from fan-shaped shells with lower angular momentum (e.g. at 8 = 30°, center). 



Equation (1) as 5 r — > 0; in other words, a piecewise function 
proportional to 1 j \Jr a — r for r < r s and zero beyond the 
caustic radius. However, we can obtain a profile consistent 
with Equation (1) in a manner analogous to the way in which 
Equation (1) was obtained, by assuming that the caustic 
is made up of many particles with slightly different caustic 
radii, normally distributed around the average r 3 . This leads 
to a Gaussian form for the distribution function in r and v r , 

/( ' r) ~16r^V2^? eXP \ 2P r 

(18)' 

A contour plot of this phase density in (r, v r ) space is shown 
in Figure 5. In Appendix C we give a derivation of this form, 
and show that integrating over all v r retrieves the functional 
form of Equation (1). 

As seen in Figure 5, this is not a perfect representa- 
tion of the phase space distribution of a caustic: the density 
contours of the model do not exactly track the distribution 
in the caustic. In particular, we obtained Equation (18) by 
assuming that the mass in the stream is evenly distributed 
along it, when in reality there is no reason to assume this is 
the case. However, both the density profile (obtained by col- 
lapsing along the v r axis) and the velocity profile (obtained 
by collapsing along the r axis) match simulated caustics 
quite well even when they are produced by mergers with 
significant angular momentum in highly nonspherical po- 
tentials. An example for ty is shown in Figure 6. There is a 
slight difference between the completely flat distribution in 
v r given by the model and the slightly peaked distribution 
seen in these two examples because the mass in the sim- 
ulated stream is not necessarily evenly distributed, but the 
difference is remarkably small. So although the details of the 
phase space distribution are not perfectly encapsulated, this 
analytically tractable model can represent the main features 
of caustics formed from mergers quite far from the ideal case 
of a radial orbit in a spherical potential. 

To keep the expression simple, we have also not taken 
into account the distribution of the material in ve and v$, 
instead assuming that all the material has zero angular mo- 
mentum so vg — = 0. In reality, even a satellite whose 
center of mass is on a perfectly radial orbit will contain ma- 



terial with a distribution of small angular momenta, and 
even in a spherical potential this will cover a distribution 
of orbital planes and tilt the distribution slightly out of the 
(r, v r ) plane. In this simplest case, the tilt results in a linear 
relationship between 9 and vg and between <j> and v^ near 
the caustic. The relationship is more complicated when an- 
gular momentum is added or the potential is not spherical, 
thanks to the precession of the orbits in these cases. Unlike 
the caustic in the r — v r plane, the particular behavior in 
these subspaces is much more dependent on the details of 
the particular interaction. In the r — v r plane the main ef- 
fect is to blur the caustic, leading to an over-estimate of the 
width S r compared to the true stream thickness. By using 
the projected r — v r space as our full phase space, therefore, 
we are mainly losing the ability to indirectly place a lower 
limit on the satellite's mass, as we discussed in Section 3. So 
this projection effect does not render too much information 
unusable. 



5 IMAGES OF CAUSTICS 

A few additional assumptions about the geometry of a shell 
will allow us to project the radial density profile given in 
Equation (1) onto the plane of the sky, allowing us in turn 
to obtain the parameters in that equation by fitting an image 
of a shell with the model. The projection process is greatly 
simplified by assuming that the material in the caustic is 
distributed evenly in angle, over some solid angle Q, s . The 
edge of the caustic can then be approximated as a spherical 
segment that spans this solid angle, and the angular extent 
of the debris can be modeled most simply as a cone. Far from 
the caustic this is not necessarily a good approximation, but 
we are only concerned with the region within a few S r of r s , 
where it appears to be adequate. In this model the shape 
and orientation of the debris in the caustic are therefore 
described by three parameters: the angles 9 3 and (f> s of the 
cone relative to the line of sight z, and the opening angle 
a of the cone (Figure 7). The solid angle enclosed by the 
cone is O s = 2n(l — cos a). The angles are defined in the 
standard way for spherical coordinates, so that (8 a ,<f) 3 ) — 
(0, 0) corresponds to a cone opening directly away from the 
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Figure 5. Contours of the model phase space density distribution for sample parameters (left); with fitted parameters overlaid on a 
caustic from simulation B; and with fitted parameters overlaid on a caustic from simulation C (right). In all panels the straight red lines 
cross at (r s ,v s ). In the center and right panels the red parabola denotes the fit used to give k, v s , and r s ; a fit to the density profile 
determined <5 r . 
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Figure 6. Probability distributions for the edges of two simulated caustics, compared to the phase space distribution integrated over 
r > r m ; n . The normalization of the model is scaled by eye to match the data; the error bars indicate Poisson error on the number of 
counts per bin. Left panel: a caustic from simulation B with r s = 90.4 kpc, S r = 0.6 kpc and r m ; n = 81 kpc. Right panel: a caustic from 
simulation C with r s = 73.9 kpc, 8 r = 1.0 kpc and r m ; n = 70 kpc. 



observer along the line of sight. If the caustic has a sharp 
edge in projection, this means that 8 S must be close to tt/2, 
in which case a is close to the value measured in the plane 
of the sky. 



For a given sky position (x, y), we then obtain the sur- 
face density E by integrating the density distribution (1) 
along the line of sight z, for z values which are inside the 



conical region occupied by the debris 
E(x,y) = A / dz 

J z m i n (6 s ,(/) s ,ot;x,y) 
-(V* 2 +!/ 2 +* 2 - 



y/ x 2 + y 2 + z 2 — r s 



x B 



Wx 2 + y 2 + z 2 - r s f 
4S 2 



(19) 



where we have defined A = fo/V^irn. The limits of the 
integration depend on the geometric parameters of the shell 
and the sky position. We can obtain an analytic expression 
for these limits, for arbitrary parameters and position, by 
means of the standard cone equations. The full derivation, 
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Figure 7. Left: the geometry of the shell and cone in the xy plane (the plane of the sky) in a slice through z = 0. The background 
shading is proportional to the density (lighter=denser), for a shell with radius r s and thickness r s /50. Superimposed are the angular 
limit of the shell a (dotted lines) and the orientation of the shell in the sky plane <j> s (dashed line). Right: the geometry of the shell 
in the xz plane (along the line of sight), sliced at y = 0. Here the background shading is proportional to the radial velocity, \v r — v s \ 
(lighter=higher value); at a given radius the radial velocity is double-valued at ±|tv — v s \. The angular limit a (dotted lines) and the 
inclination B s with respect to the line of sight (dashed line) are superimposed. A line of sight (red solid line) with a given projected 
radius R, equal to x in this slice, probes a range of radial velocities and therefore line-of-sight velocities v z . The range of v z is limited 
by the values of v s and ft according to Equation (23), but also by the angular extent of the shell and its orientation with respect to the 
observer. 



and the method for calculating z m i n and 2: max , are given 
in Appendix D. The integral along the line of sight must 
be done numerically, but since the expressions for z m i n and 
Zmax are analytic in 9 S , cj> 3 , and a, standard minimization 
routines can be used to find best fit values by comparing 
the calculated and measured profiles. To this end we also 
provide derivatives of the integration limits with respect to 
the parameters in Appendix D. 

Two examples of the result of this process are shown in 
Figure 8. We used the N-body realizations of two caustics 
from simulation D (left column) to construct two binned 
surface-density distributions (center column), which were 
then least-squares fit to Equation (19) to obtain the pa- 
rameters (r 3 ,8 r ,A), which describe the radial density dis- 
tribution, and the parameters (0 3 , <f) 3 ,a), which describe the 
orientation and angular extent of the shell. The best-fit mod- 
els are shown in the right-hand panel. Although the data 
clearly have smaller-scale complexity that is not accounted 
for in the model, the model is approximately consistent with 
the smooth features of the shell when comparing by eye. 
The surface-brightness fit is mainly used to provide an in- 
put geometry for fitting the kinematic data, from which g s 
is recovered, so a model that roughly describes the extent, 
variation, and orientation of the shell is sufficient. 

One significant exception is seen in the top row of Fig- 
ure 8. Compared to the model, the caustic taken from the 
N-body simulation appears to have r s smaller toward the 
top of the frame. Correlations between the energy and or- 
bital phase of the material in the stream, combined with the 
presence of angular momentum, will introduce a variation in 



the caustic radius with physical angle r s (#, <j>) which is not 
allowed for in the model. This variation can be comparable 
to S r . This discrepancy can be addressed by allowing r s to 
vary, but since 5 r is mainly used as input into the kinematic 
models, this extra parameterization may not be needed. 

The normalization A involves the phase-space density 
/o and the gravitational force at the caustic g B . If the model 
is fit to a surface- brightness map, the normalization obtained 
by the fit must be scaled by the estimated mass-to-light 
ratio of the debris to obtain the normalization of the mass- 
density profile. In principle, if we use Equation (11) as an 
approximation for k and guess a reasonable value (or range) 
of /o, we can obtain an estimate of g„: 



9a 



nX 2 
" fS 



(20) 



However, due to the strong inverse scaling of g s with /o, the 
estimate of the phase-space density has to be within better 
than a factor of about 3 just to get the right order of magni- 
tude for the gravitational force. Furthermore fo refers to the 
fine-grained density, not the coarse-grained density probed 
by the stars, so it is unlikely that this approach will yield a 
reasonable constraint on g s . We conclude that realistically, 
images alone can only constrain the combination fo^/gl. As 
we discuss in the next two sections kinematic data can sep- 
arate this constraint into independent measurements. 

Although it cannot produce an independent measure- 
ment of g s , the surface-brightness fit plays a crucial role 
in the interpretation of kinematic data. Fitting the image 
models the geometry of the caustic, determining all the pa- 
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Figure 8. The binned "images" (relative brightness) generated 
for views of two caustics from simulation D (left: points from 
Nbody simulation; center: binned surface density) can be ade- 
quately modeled by a projection of Equation (1) onto the sky 
plane using a spherical cone to model the geometry as described 
in the text (right). The caustic parameters (r s , <5 r , A., a, 9 3 , rf> s ) 
used to generate the model image on the right in each case were 
obtained by a least-squares fit to the mock data in the center 
image. 



rameters of the phase space distribution except for k and 
M s . These two parameters are necessarily convolved with the 
shell geometry in line-of-sight velocity measurements, so de- 
termining the geometry independently of the kinematics al- 
lows the most efficient use of the spectroscopic information 
which is so difficult to obtain. In the next two sections we 
discuss two ways that line-of-sight velocity data could be 
used to estimate g s . 



6 DISCRETE LINE-OF-SIGHT VELOCITY 
MEASUREMENTS OF CAUSTICS 

For some shells, the line-of-sight velocities of point sources 
in the shell can be measured. These measurements probe 
the projection of the phase space distribution of Equation 
(18) onto the space of projected galactocentric distance R 
and line-of-sight velocity v z . As was shown in Merrifield & 
Kuijken (1998), this projected distribution can have a dis- 
tinctive, roughly triangular shape for shells that are oriented 
nearly edge-on with respect to our line of sight. The shape 
of the outer envelope of the profile, like the surface bright- 
ness profile, depends on a few parameters that can be fitted 
to obtain constraints on the gravitational field at r s . More 
specifically, these parameters determine the relationship be- 
tween the projected radius R and the maximum magnitude 
of the line-of-sight velocity at that radius, w z , max . Here we 
derive an expression for u z . max (-R), similar to the derivation 
in Merrifield & Kuijken (1998) but properly including the 
nonzero apparent outward velocity of the caustic. We then 
compare the result to simulations, and test how well the fit 
of this function can recover g B . 



6.1 Derivation of the envelope equation 

A given projected radius R samples different three- 
dimensional radii r that correspond, via the phase space dis- 



tribution, to different line-of-sight velocities v z . For a satel- 
lite on a purely radial orbit, 



Vz =V r -, 

r 



(21) 



so that for a given r and R, the line-of-sight velocity at the 
midline of the stream is described by 



v s ± 



vV - R 2 



(22) 



To find the maximum v z at a given projected radius R we 
need to maximize this expression with respect to r, to find 
the physical radius r e that contributes the highest line-of- 
sight velocities. Technically this results in a 6th-order poly- 
nomial expression for r e ; however the approximation used 
by Merrifield & Kuijken, r e w (r s + R)/2, works quite well 
for the region near the shell. 

Using the linear approximation for r e and the first-order 
approximation for « gives the formula for the maximum line- 
of-sight velocity as a function of projected radius alone, with 
g(r s ) and v s as parameters: 



vVf + 2Rr s - 3i? 2 (^/g{r s ){r s -R) + i>„) 



+ R 



(23) 



Taking the limit v s — > and keeping only first-order terms 
in an expansion in R about r s , we can recover the expression 
obtained by Merrifield & Kuijken: 



R). 



(24) 



We can derive an expression for the density in the pro- 
jected phase plane, f(R,v z ), by converting the full / to 
cylindrical coordinates and integrating over everything but 
R and v z to project the distribution onto the plane of the 
sky and the line of sight. This means we must take the ge- 
ometry of the shell into account, using the cone model de- 
rived in the previous section. Schematically, the full phase 
space distribution is limited to the cone by some function 
0(9 s ,(f> s ,a;9,(f>) in spherical coordinates. Since we are pro- 
jecting onto the plane of the sky, we can separate the in- 
dependence so that the cone limits are 



e(e s ,4> s ,a;R,z)@(\(f>-(t> s \ < a). 



(25) 



We have also assumed the velocity of the material is purely 
radial, and so ve and v$ are both zero. In cylindrical coor- 
dinates this becomes 



5(v^, — 0)S(vr = v z R/z). 
Then the full phase-space distribution is 



(26) 



f(x,p) = .Foexp- 



VR 2 



"R _|_ 



xB(O a ,<l> a ,a;R,z)e(\<j>- <£,| < a) 
x<5(« = Q)S(vr = v z R/z) 



where 



To 



15 



(28) 



327rr? /2 (l - cos a) V 27tS ? 
Now we must integrate over v^jURjCj), and z. All of these 



(27) 
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Parameter 



Value, Caustic 1 Value, Caustic 2 



r s (kpc) 92.8 

5 r (kpc) 0.32 

v s (km s- 1 ) 67. 

| Ss | (km s- 1 Myr- 1 ) 1.2 



47.7 
0.15 
40. 
4.0 



Table 2. Parameters for the phase-space models of the two caus- 
tics shown in Figure 9, derived by fitting the full phase-space 
information as described in the text. 



but the z integral are delta functions or step functions, and 
can be trivially evaluated: 



This is a result of correlations in the perpendicular 9 — v$ 
space that are specific to spherical potentials. The velocity 
vg varies linearly with 6 along the stream thanks to the sym- 
metry in this coordinate, is centered on zero near the caustic, 
and has a range of values comparable to v s in the range of 
interest in R. For this highly symmetric system, this corre- 
lation causes a rounding of the outer envelope and leads to 
a slight overestimate of g. If we examine a case where the 
angular momentum is nonzero (Simulation B, Figure 10), 
we see that the addition of angular momentum breaks the 
symmetry in the 6 — vg plane, leading, counterintuitively, to 
an improvement in the fit. 



f(R,v z ) = 2T a / dz<d(6 s ,<j> s ,a;R,z) x 



exp • 



[r s - VR 2 + z 2 - + v s ) 2 ]' 

2S 2 



The remaining step function limits the line-of-sight integral 
to points inside the spherical cone used to describe the shell's 
geometry; thus we can rewrite the integral by replacing the 
limits of integration with the functions of R derived in Ap- 
pendix D (suppressing the parameter dependence for sim- 
plicity), just as we did to calculate surface brightnesses: 



/•Zmax(fl) 

f(R,v z ) = 2J r oce / dze 



The line-of-sight integral is then performed numerically to 
obtain the distribution for a given R and v z . 



6.2 Comparisons with simulations 

To see how well the envelope equation and the projected dis- 
tribution in R and v z represent the main features of shells, 
we compare them to simulated tidal caustics. We first con- 
sider simulation A, which satisfies all the assumptions used 
to derive the models. The top panel of Figure 9 shows the 
system in the (r, v r ) projection of phase space, with two 
caustics identified in red and blue. The simulation provides 
us access to the full phase-space information of the shell: as 
a first test of the model's ability to represent the observed 
quantities, we can get the parameters of «z, m ax(-R) from the 
complete phase space information and compare them in the 
observed space (R,v z ). To determine the parameters we fit 
Equation (6) to each caustic (solid cyan and magenta lines) 
to obtain values for r s , v s , and g s ; we also obtain 8 r by fit- 
ting the radial density profile of each caustic (a histogram 
of the r coordinates of the colored particles) with Equation 
(1). We need 5 r because the spread around r s is not taken 
into account by Equation (23), so the envelope should be 
shifted radially outward to the very edge of the distribution 
at r s + 28 r . The parameters are given in Table 2. 

The center panel of Figure 9 shows the same system in 
the observed space (R,v z ). For both caustics it is clear that 
Equation (23) (cyan/ magenta) is a better description of the 
outer envelope than Equation (24) (green). However, some 
material still lies above the envelope described by the equa- 
tion for Wmax, especially for the outer caustic, even though 
we have adjusted for the spread of the material around r s . 



The lowest row of Figure 9 shows the modeled projected 
density in the R — v z plane, Equation (30), for each of the 
2 "jtwo modeled caustics. The model captures the interior struc- 
ture in the density distribution that is seen in the N-body 
[distribution (second row), including significant density en- 
chancement in the interior of the outer (red) caustic in the 
(2^rojected space. In practice, the number of measurements in 
this space is likely to be fairly small, probably insufficient to 
distinguish the structure of the distribution by eye. Because 
the outer envelope is not filled with a constant density of 
points, a sparse sampling of measurements in this plane will 
not give a good estimate of t) z , max (fl). Fitting Equation (23) 
directly to an ensemble of measured line-of-sight velocities 
is therefore not a good way to determine the gravitational 
j^r s - v / 'r, 2 +z 2 -k( r ^ z +^-u a ) 2 j /2iforce. Instead, the measurements should be assigned proba- 
bilities for a given set of model parameters using Equation 
(30). Then the parameters, including g, can be recovered by 
maximizing the total likelihood. 



(30) 



Next we consider an example in which the potential has 
a realistic axisymmetric contribution from a disk: simulation 
D (Figure 11). (We will return to simulation C later in this 
section.) The caustics in simulation D combine three im- 
portant differences from the highly symmetric situation in 
simulation A: they are not aligned exactly with the line of 
sight, they were formed with some small angular momen- 
tum, and the potential is not spherical. In the plane of the 
sky (leftmost panel), the outermost caustic, shown in blue, 
appears less sharp-edged than the next caustic, shown in 
green. As before, we fit the r — v r profiles of the caustics 
(middle panel) to obtain r s and v s ; in this case we also use 
the fitted value for k which is not far from the true value of 
g(r s ). In the projected space (right panel) we see again that 
the velocity profiles are not all symmetric around v z = 0, 
but the displacement can be adequately modeled by a con- 
stant, and the distributions are fairly symmetric around an 
offset velocity axis. We will show momentarily that these are 
primarily projection effects. Second, the outermost caustic is 
rotated enough to affect the velocity profile, since the bulk 
of the material does not satisfy the approximation R « r 
(or 6 S « 7r/2) used in the derivation. Nonetheless, the enve- 
lope fits the material that does satisfy this approximation, 
and when the envelope is translated in R to encompass the 
bulk of the observed material (dashed black line) it is still 
a roughly accurate description. The difference between the 
very outer edge of the shell and the outer edge of the main 
part of the material is about 5 kpc, so the uncertainty in 
R will not result in a large error in g s . This suggests that 
Equation (23), with some modifications, should be able to 
recover g s for real shells. 
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Figure 9. Line-of-sight velocity profiles and phase space distributions for two caustics from simulation A. The i — v r projection and the 
selected caustics are shown in the top panel (blue and red points) with quadratic fits overplotted (cyan and pink lines). The middle panel 
shows the projected line-of-sight velocities in the plane of the sky with Equations (24) (green) and (23) (cyan and pink) overplotted, 
using the values for r s and v s obtained from the fits in the top panel and k = l/|2</(r s )| for the input potential. The bottom panel shows 
the same two equations (green and blue/red dashed, respectively) superposed on the projected density, Equation (30), for the same two 
caustics. 
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Figure 10. As in Figure 9, but using simulation B. The nonzero angular momentum actually improves the agreement of the model and 
simulation in the (R,v z ) plane somewhat. 
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Figure 11. Views similar to Figure 9, but for simulation D. The red, green, and blue regions are the same as in Figure 1. The left panel 
shows the view in the r — vr plane with fits to the colored regions as black lines. The right panel shows the R — v z plane with envelopes 
calculated using Equation (23) and the values from the fits in the center panel; the two inner envelopes are translated by eye in the v z 
direction to fit the symmetry axis. The dashed line is a translation in R of the envelope for the outermost caustic. See the text for further 
discussion. 
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6.3 Projection effects 

The orientation of the shell with respect to the line of sight 
can affect the view in the (R,v z ) plane. In Figure 12 we 
show the effect of rotating the outermost caustics from sim- 
ulations A and B from a perfectly symmetric orientation 
(6 3 = 90°) away from the observer (6 S < 90°). All these ex- 
amples have 6 3 very different from 90°; at inclinations closer 
to 90° there is not very much change to the distribution, but 
if the shells are nearly spherical (as these are) then even at 
relatively large inclinations they may still have a sharp edge 
in projection. We see that there are three main effects: the 
peak density is no longer on the v z axis, the distribution 
is no longer symmetric about the v z value defined by this 
point, and the internal structure seen in the distribution is 
smoothed out. 

The distribution is offset in v z because the average dot 
product with the shell expansion velocity v s , and any net 
rotation in the line-of-sight direction, is now nonzero. In the 
figure we have modeled this shift by 

«offsot = (v s + cosd s , (31) 

where the first term accounts for the shift from the shell 
velocity and the second for the rotation. As shown in the 
figure by the dashed lines, this model is a good predictor 
of the vertical shift for spherically symmetric potentials. In 
the case of simulation B, the angular momentum is in the 
—z direction so that the two terms nearly cancel. This also 
causes the lower half of the distribution to be selected in- 
stead of the upper half for this case, even though in both 
simulations the viewing angles are the same. 

The asymmetry about the v z — foffset axis is a combined 
effect of the inclination angle and the finite angular extent of 
the shell — essentially, the line of sight now probes an asym- 
metric distribution of velocities. However, it is clear from the 
figure that the "filled" side of the distribution still has the 
same envelope as the symmetric distribution. Rotating the 
caustic in the other direction (toward the observer instead 
of away from them, 6 S > 90°) will "fill" the other side of the 
distribution instead, but the maximum value of \v z — t) ff ae t| 
still follows the model. Additionally, for caustics inclined to 
the line of sight the projection smoothes out the interior 
structure of the distribution somewhat, which could be an 
advantage when sampling this plane with discrete sources 
since the sampling is then more uniform in the populated 
part. In any case, the only asymmetry that the projection 
effect can produce is to cut out part of the filled area under 
the envelope. This may complicate efforts to fit the enve- 
lope, so an understanding of the geometry, especially the 
inclination angle 6 3 , is important since kinematic measure- 
ments are unlikely to populate the whole (R,v z ) plane the 
way these simulations do. 

6.4 Indications of a flattened potential 

A significantly flattened potential like that of simulation C 
can complicate the shape of the (R,v z ) plane in ways that 
cannot be obtained by inclining a shell to the line of sight 
as discussed in the previous section. For extremely flattened 
potentials the kinematic plane can look quite different from 
the examples presented above, even if the caustics appear 



shell-like in position space as they do in this case (see Figure 
1). This example is intended to be especially tricky because 
the potential is not only extremely flattened, but face-on 
(the symmetry plane is the x-y plane) so that the degree 
of flattening could not be estimated from the galactic shape 
for this case. 

The caustics from simulation C, shown in Figure 13, 
do form the characteristic structure (though significantly 
broader) in the (r,v r ) plane (top left), but the (R,v z ) pro- 
jection (top right) now looks remarkably different. The caus- 
tics are not symmetric about v z = 0, as expected for shells 
inclined relative to the line of sight, but additionally the 
slope of w z ,max(i?) no longer agrees with the model predic- 
tion based on a fit to (r, v r ). The solid lines in the top right 
panel of Figure 13 show Equation (23) shifted by eye in v z 
to the rough axis of symmetry of each shell, and show that 
this approximation is not a good description of the distri- 
bution in the projected phase space. This failure is due to a 
combination of effects. First, examining the phase space in 
projections consistent with the potential's symmetries (Fig- 
ure 13, bottom row) shows that especially for the outer (red) 
caustic, the material is actually passing through two caus- 
tics, one in R and one one in z. So whereas in a spherical 
potential v z is single-valued with respect to z, thanks to 
this double caustic structure it is now double- valued. Fur- 
thermore, the caustic in (R, vr) is by definition not probed 
by the line-of-sight velocity, where the spherical caustics had 
a component of v r along the line of sight. In addition, we see 
from the range of z values in the lower right panel that the 
assumption that R « r is no longer valid. As a result the 
dot product of the velocity with the line-of-sight direction is 
now primarily giving us information about the behavior of 
v z , not the apparent caustic in (r,v r ). 

We include this example to demonstrate how measure- 
ments of the kinematics in caustics can diagnose a signif- 
icant amount of flattening in the potential. The degree of 
displacement of the (R,v z ) distribution from the v z axis 
shown in Figure 13 is quite large compared to those ob- 
tained by rotation in Figure 12; interpreted in the context 
of this model, one would have ventured that the distribution 
seems to indicate a large amount of angular momentum com- 
bined with a 6 S very different from n/2. However, the shells 
still look sharp in projection on the sky, which is inconsis- 
tent with a large inclination and unlikely for such a large 
L. The distribution in (R,v z ) very clearly does not follow 
the model envelope on either the positive or negative side, a 
type of asymmetry which cannot be replicated by adjusting 
the viewing angle, and includes extra lumps and curvature 
that the model does not. These signatures, whether probed 
by discrete measurements or integrated spectra, could thus 
potentially be used to diagnose departures of the potential 
from spherical symmetry. 

6.5 Recovery of the gravitational force 

We can also use the simulations presented here to quanti- 
tatively test the ability of this method to recover the value 
of g s in the case of projected observables. For each caustic 
we need to fit the envelope of the points from the N-body 
simulation in the (R, v z ) plane. To do so, we bin the points 
in R and choose « Zjma x in each bin at the nth percentile in 
v z , where n is usually between 95 and 100. We choose n as 
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Figure 12. Effect of shell orientation relative to the line of sight on the (R, v z ) plane, in simulations A (left) and B (right). The colors 
red through purple indicate 8 a = {90,63,60,53,49,45}°. The dashed black lines are the envelope derived in the text, translated in v z 
by the amount defined in Equation (31). The v z shift and asymmetry are results of taking 8 S 5^ 90°, but the envelope still fits the more 
fully populated side of the distribution. 




Figure 14. Example of envelope data and fits for the outer caus- 
tic from simulation B. The black points with error bars are ob- 
tained from the red points (the simulation particles) as described 
in the text; the green line is the fit to Equation (24) and the 
pink line is the fit to Equation (23). This is one of the less well- 
fit cases, as indicated by the overshoot on the front edge of the 
distribution, which leads to a larger uncertainty on the recovered 
9s- 



large as possible in each case so that we get the edge of the 
envelope but are not contaminated by single outlier points 



above it. We estimate the errors on these v z 



points by 



looking at the difference in v z between roughly the n and 
n — 2 percentile, to get a sense of the "fuzziness" of the edge 
in each bin. Finally, we add a point at the front edge of 
the distribution where u z , ma x goes to zero. This method is 
not meant to be representative of how real data would be 
analyzed (for one thing, it is not likely that we would have 
thousands of real velocity measurements to work with) but 
was developed as a rough way to represent the envelopes 
of the particle representations of the simulated caustics. An 
example is shown in Figure 14. We then carry out weighted 
least-squares fits of Equations (23) and (24) to recover g s , to 
test whether one form is better than another at measuring 
g s and how close we can get to the input value. 

The results, shown in Figure 15, show that fitting Equa- 
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Figure 15. Recovered values of g s for the caustics from simula- 
tions A-D pictured in Figures 9-11 (black=A, red=B, green=C, 
blue=D). Each caustic was fit twice, once using Equation (24) 
(points with dashed error bars) and once using Equation (23) 
(filled squares with solid error bars). The error bars are the pro- 
jected 95 percent confidence intervals from the least-squares fits. 



tion (24) (open squares) consistently overestimates g s by 
large factors, up to a factor of 4 in some cases. Fits to Equa- 
tion (23) (filled circles) are consistently better at recovering 
g s , although in a few cases the enormous error bars indi- 
cate that this can be a more difficult fit to converge. The 
only cases in which both equations are more or less equally 
bad are the two caustics from simulation C (green points), 
which we do not expect to give a good result for the reasons 
discussed above. In most cases we can recover g s within 50 
percent using Equation (23). 



7 SPECTROSCOPIC LINE-OF-SIGHT 
VELOCITY MEASUREMENTS OF 
CAUSTICS 

Spectroscopic measurements of integrated light use a slit 
or a set of fibers to record the line-of-sight velocity profile 
at a given location in the plane of the sky. We can calcu- 
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Figure 13. Top row: as in Figures 9 and 10, but for simulation C. Here the model, which assumes that the shell lies around the plane 
of the sky so that z is small, does not work even when translated in v z to the axis of symmetry for each caustic. Bottom row: slices of 
the phase space in natural axisymmetric coordinates show that material is actually passing through two fold caustics at once, in Ft (left) 
and z (right). Also note from the z axis that the assumption that r ~ R used to derive Equation (23) is not valid for this caustic. 



late the expected velocity distribution in a slit or fiber from 
the phase-space distribution. As an example, we consider 
a single round fiber of the type that might be used in an 
integral-field unit. The center of the fiber is located at some 
sky location (x&b,yRb) and has a radius pah- We take these 
quantities to be in physical units (e.g. kpc) so that g(r s ), 
embedded in the phase space distribution, still has reason- 
able units. To get angular units for the fiber values one has 
only to divide by D, the distance from the observer to the 
target. The supplied profiles can then be convolved using 
spectral simulator tools to produce the expected lineshape 
for fitting purposes. 

For a circular fiber of arbitrary size, the velocity profile 
requires three integrals: over the line of sight, and over aux- 



iliary variables £ and r\ that span the area subtended by the 
fiber: 



f(v x ) = 2F a dU dr) dz 



exp 



where 



and 



R= \A%sb + £) 2 + (ytib + V) 2 



(33) 



r = \/{xRb+0 2 + (VRb + v) 2 + z 2 - (34) 
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For a different spectroscopic geometry, e.g. a slit instead of 
a fiber, one would simply alter the integration ranges of £ 
and rj in this expression. 

To simplify the calculation we can assume that the an- 
gular size of the fiber is small compared to the angular size 
of the shell, so that it probes a range of physical R compara- 
ble to or smaller than the shell thickness S r . For many shells 
this is likely to be a sufficient approximation; for example, 
a 4" fiber has a fiber radius of about 0.5 kpc at a distance 
of around 25 Mpc, comparable to 5 r for the shells in our 
simulations (see Table 2). Assuming that each fiber probes 
a single line of sight allows us to replace the limits of the 
line-of-sight integral with their central values z^™ (xfib, 2/flb), 
and the £ and rj integrals with 7rp| b times the central value 
of the integrand. This reduces Equation (32) to a single in- 
tegral over the line of sight: 



f{v z ) = 2irp fih J r a J i/-. cxp<J- 
R& b v z 



+ 



(35) 



where 



Rtib 



+ vL- 



(36) 



For some lines of sight this distribution has four sym- 
metric peaks, while others closer to the shell have a double 
peak (Figure 16). The shape of the profile depends on the 
distance from the shell edge, the shell thickness, the orienta- 
tion of the shell relative to the line of sight, and the angular 
span of the shell. The distance from the shell edge deter- 
mines whether two peaks or four are observed: for the line 
of sight in Figure 7, ignoring the angular cutoffs momen- 
tarily, four peaks will be observed because the line of sight 
probes radii far from the shell edge, where the two values of 
±|« r — v s \ are widely separated, closest to z — 0. The line of 
sight probes both near and far edges of the shell, giving rise 
to two mirror-image pairs of peaks. For a line of sight closer 
to r s , however, v r is closer to v s , so the two velocity val- 
ues at a given radius are close together and each bifurcated 
pair of peaks will merge into a single one. The distance from 
the shell edge where this merging happens depends on the 
shell thickness, but also on fiber size and spectral resolution. 
These spectral line shapes are similar to those presented in 
Ebrova et al. (2012), but include the broadening and changes 
in peak position induced by the range of energies of the stars 
in the shell. Their velocity profiles can be obtained from our 
expression by taking the limit 5 r — > 0, or equivalently us- 
ing Equation (17) for the distribution function instead of 
Equation (18). 

The locations of the centers of the outer peaks are given 
by (R) as in Equation (23), but since the velocity pro- 

file convolves the physical density of the material with its 
shape in (r,v r ) space, we should replace r s in the equation 
with r max , the radius of peak density defined in Equation (3). 
Thanks to the stream thickness these two radii are slightly 
offset from one another. The outer peak locations computed 
this way are marked in the left-hand panel of Figure 16; 
the colors match the respective lines of sight marked in the 
center panel. The interior peaks' locations can likewise be 
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Figure 18. Changing the value of g(r s ), for an otherwise con- 
stant set of shell parameters and constant line of sight (along 
the symmetry axis of the shell near its edge), alters the veloc- 
ity profile. The colors from red to purple are the values g(r s ) = 
{1.2,1.5,2.5,3.4,4.9,7.6} km s" 1 Myr -1 , corresponding to the 
gravitational force at r = {90, 80, 60, 50, 40, 30} kpc in the po- 
tential of simulations A-C and the range r = 10 — 35 kpc in the 
potential of simulation D. Stronger gravity results in a broader 
profile and lower peaks. The shell is the same model as in Figure 
16; the orange profiles in the two figures correspond. 



computed by letting v s — > —v 3 in Equation (23) and again 
using r 3 — > r max . So one way to measure g s from spectra is 
similar to that discussed in Section 6 for discrete measure- 
ments, except that since the spectra probe the integrated 
light they will automatically measure the correct v max - How- 
ever, the caveats about potential flattening and angular span 
discussed above still apply: for this reason it is important to 
understand the geometry of the shell by fitting the surface 
brightness profile in order to judge whether Equation (23) 
is applicable. 

The inclination and angular span of the shell have an 
effect on the integrated velocity profile just as they do on 
the distribution in (R,v z ). Considering the angular cutoffs 
in Figure 7, we see that the limited angular span of visible 
material restricts the velocity profile so that only one of 
the two pairs of peaks (one of the two intersections with 
r s ) is probed. For S near tt/2, the line of sight is roughly 
symmetric around z = so the two pairs of peaks are both 
the same height, but for inclinations much different from 
this, lines of sight can intersect only some of the peaks and 
cut off others, as shown in Figure 17. Knowing the geometry 
of the shell beforehand from image fitting is therefore a great 
advantage when deciding where to place fibers, since the 
chosen lines of sight must intersect all the peaks in order to 
measure v x , max . 

The profile in even a single fiber is sensitive to g(r s ) 
(Figure 18), so g s can in principle be recovered by fitting 
the full velocity profiles of even a handful of lines of sight, 
even if there are too few to get the variation of the peak 
widths and fit Equation (23). Fitting the lineshape also al- 
lows the combination of multiple fibers (for example, from 
an integral-field unit) to improve the data quality. We will 
assess the quantitative ability of this strategy to recover g s 
in future work. 
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Figure 16. Equation (32) can be implemented numerically to predict the velocity distribution for different lines of sight through a shell. 
Here we illustrate using a particle realization (left) of a shell modeled on the one shown in green in Figure 11. The different lines of sight 
(center; colored circles) are overlaid on a surface-brightness map of the model convolved with a Gaussian of 2r^; the red box in the 
right-hand panel indicates the area shown in the center panel. This example uses a fiber radius of 0.5 kpc, equivalent to 4 arcsec at 25 
Mpc, shown to scale in the colored circles of the center panel. The different lines of sight produce different velocity profiles (right) based 
on their distance from the shell edge. Farther from the edge the profile has two pairs of peaks, while close to the shell edge only one pair 
of peaks is visible. The colored vertical lines mark the outer peak locations predicted by Equation (23) as described in the text. 
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Figure 17. The precise angle of the shell with respect to the line of sight influences the velocity profile. In this example, as the shell is 
inclined more and more away from the observer the line of sight (black dashed line in the left panel) starts to probe a smaller range of 
radii. In this panel the shell occupies the intersection volume between the cone and the sphere. The angular cutoff causes an alteration 
in the velocity profile as a smaller range of velocities is included (right). The right panel shows the integrated velocity profile along the 
same line of sight for the inclination angles 8 S = {90, 63, 60, 53, 49, 45}° for red through purple; the colors in the two panels correspond 
with each other and with the colors in Figure 12. 



8 DISCUSSION AND CONCLUSIONS 

In this paper we have presented a simple analytical model 
that successfully describes tidal shells with a handful of pa- 
rameters, some of which are linked to the masses of the two 
interacting galaxies. The model includes self-consistent ex- 
pressions for the surface brightness (Sections 3 and 5) and 
phase-space distribution (Section 4) of material near the 
edge of a single tidal shell. Given an image of a shell, the 
surface-brightness model jointly constrains the fine-grained 
phase-space density of the satellite material, /o, and the 
gravitational force at the shell edge through the parame- 



ter k, which is related to g s as discussed in Section 3. The 
gravitational force g s (r s ) can be independently measured if 
line-of-sight velocities are also obtained. These can be either 
the velocities of individual point sources in a shell (Section 
6) or velocity profiles obtained from integrated-light spectra 
(Section 7). If multiple shells are observed around a galaxy, 
each one gives an independent measurement of the gravita- 
tional force at the radius of its edge. 

The model makes three assumptions about the two in- 
teracting galaxies: that the interaction orbit is radial, that 
the potential is spherical at the shell radius with negligible 
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tides, and that the satellite originally had a Maxwellian ve- 
locity distribution. We demonstrate in Section 3 that the 
assumptions of zero angular momentum, zero tides, and 
spherical symmetry in the potential are sufficient for caus- 
tics that appear as fan-like shells on the sky. The model is 
self-consistent in its treatment of the effect of the spread of 
orbital energies in a given shell, except for the assumption 
that the shell radius is constant. The model does not allow 
for a non-uniform distribution of material in the stream. 
In Section 4 we show that variations in the bulk velocity 
profile from the slightly nonuniform density along simulated 
streams are minor. 

From the analytical model, we calculate different ob- 
servable quantities — images, line-of-sight velocities, and ve- 
locity profiles — and demonstrate how each observable is sen- 
sitive to the radial component of the gravitational force, 
g s , acting on the debris at the caustic. Thanks to the rela- 
tively simple forms of the phase space and spatial density, 
the calculations can be easily tailored to a specific observ- 
ing paradigm. In Section 5 we show an example of how a 
simple geometric model, combined with the density profile, 
can be used to construct a basic but sufficient representation 
of simulated shells. Although images alone can only jointly 
constrain g B and /o, the image fit is important to determine 
the geometry of the caustic, which is imprinted on line-of- 
sight velocity measurements. 

In Section 6 we show that the maximum-velocity en- 
velope of pointwise line-of-sight velocity measurements can 
be used to recover g 3 to within better than 50 percent in 
most non-pathological cases; the crucial ingredient is to al- 
low for the nonzero expansion velocity at the shell, which 
is a reflection of the energy spread of material in the caus- 
tic. The kinematic model fails when the potential is so flat- 
tened that it produces caustics in more than one dimension 
nearly simultaneously; this can be diagnosed by examining 
the symmetry and substructure of the velocity distribution 
around the mean v z . The viewing angle of a shell can also 
produce an asymmetric velocity distribution, but only in a 
specific way; it cannot mimic the effect of a significantly flat- 
tened potential. The peaks of the velocity profiles presented 
in Section 7 depend on g a in the same way as the velocity 
envelope, with the added advantage of probing the entire 
distribution along a given line of sight rather than sampling 
individual points. The value of g s can also affect individual 
lineshapes, which could lead to another way of measuring 
this quantity. 

One possible caveat to this work is that even for high- 
mass-ratio mergers like those studied here, the use of a static 
potential is not necessarily justified during the disruption of 
the satellite at pericenter, where the enclosed mass of the 
host galaxy can be comparable to the satellite mass. We con- 
ducted a high-resolution simulation using a live mass distri- 
bution, with conditions similar to simulation B, to study the 
effect of the response of the inner host galaxy. In agreement 
with previous work (e.g., Seguin et al. 1996), adding this 
response results in slightly thicker, more massive caustics at 
slightly smaller radii, but does not destroy the shells' phase 
space structure. In addition, Mori & Rich (2008) performed 
Simulation D with a live representation of M31's potential 
and found that the same shells-and-stream system formed in 
that case. This suggests that shells are not fundamentally 



different in live host galaxies and that our results should 
extend to these cases. 

Finally, each shell observed around a galaxy provides 
its own mass estimate at its own radius, without having to 
assume anything about the relationship between different 
shells. In galaxies with many shells, therefore, this technique 
could be used to obtain a rough dynamical mass profile of the 
host galaxy at otherwise inaccessible locations outside the 
luminous disk or ellipsoid. By using their distinctive mor- 
phology as a clue to their symmetry, shells will open for us 
a new window into the shapes and masses of external galax- 
ies. 
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APPENDIX A: POTENTIALS USED IN 
SIMULATIONS 

Simulation A and Series B (see Table 1) use the spherical 
isochrone potential with total mass M and scale radius b 
described by 

GM 



$ sph (r) = -- 

iso V / 



(Al) 



b + Vr 2 + b 2 ' 

where r is the spherical radial coordinate and G the gravi- 
tational constant. 

Series C uses a flattened axisymmetric version of the 
same potential, 

GM 



b + sjB? + z^/q 2 + & 2 ' 



(A2) 



where q is the flattening parameter. In the limit q — > 1 
one recovers the spherical isochrone potential. For simula- 
tion C we take q = 0.7 which produces an extremely flat- 
tened distribution; contours of the potential and radial force 
are shown in Figure A. 

Simulation D uses the combined potential described in 
Table 2 and the references therein. 



APPENDIX B: DERIVATION OF THE 
NORMALIZED COLD PHASE SPACE DENSITY 

Near a caustic, the phase-space distribution function (DF) 
can be approximated by a line, a quadratic equation relating 
the galactocentric radius r and the radial velocity v r : 

\21 



f(r, v r ) = FqS [r s - r - k (v r - v s ) 2 ] 



(Bl) 



Here we have suppressed notation indicating that there is no 
variation in ve or (two trivial delta functions) and also 
notation limiting the angular extent of the shell to some to- 
tal solid angle f2 s of the sphere. We will deal with this second 
notion by adjusting the limits of integration in angle. To de- 
termine the normalization To we impose the normalization 
condition 



fr s f fo 

/ r 2 dr dQ 
Jo Jn s J~< 



dv r dvgdv$f{r, v r 



(B2) 



The integrals over ve and v$ cancel, and the integral over 
the solid angle spanned by the shell, fl s , simply contributes 
a prefactor. Thus this expression simplifies to 



(B3) 



rr s poo 

FqVLs I r 2 dr / dv r S [r 3 — r — n (v r — v 3 ) 2 ] = 1. 

JO J-oo 

We now change variables to the argument of the delta func- 
tion, 



K (v r — V s ) = A r — K (v r — V s ) 



(B4) 



to integrate over v r . Note that the quantity A r = r s — r is 
always positive since the caustic only includes material at 
r < r s . 



This change of variables has two solutions for v r 

I A r — u 

Vr.± = V 3 ± 



(B5) 



To understand how to do the integral we therefore make a 
preliminary auxiliary change to the variable x = v r — v s , so 
that the integrand is symmetric in x: 



/oo 
dxS(A r 
-oo 



In terms of u, 



x± = ± 



and u = A r — nx 



(B6) 



(B7) 



Now we break the integral up into positive and negative x 
pieces, so that it is clear which branch of the solution to take 
on each piece: 

/{) poo 
dx S(A r — kx 2 _) + / dx S(A r — kx\) (B8) 
-oo Jo 

We note from the change of variables that at x = ±oo, u = 
— oo, while at x = 0, u takes on its maximum (positive) value 
of A r . So the zero point will still be included in both pieces 
of the integral, allowing the delta functions to be evaluated. 
To complete the change of variables we need the Jacobian 
dx/du: 



dx± 
du 



1 



(B9) 



so that 



- - r 

J — oo 

- f 

J — oo 



2yJti(A r - U) 

du 

\A(A r - 



2^/k(A t - u) 

5{u) - [ °° - du S(u 

JA r 2y/K{A r — U) 



8{u). 



(B10) 



Now we can evaluate the integral directly thanks to the delta 
function: 

K = . = = ■ (Bll) 
^Jn(r s - r) 

Note that this integral gives the functional form of the den- 
sity distribution, which has the characteristic l/y/r falloff 
behind the caustic. This shows that we have chosen an ap- 
propriate form for the DF, since it gives the right density 
distribution when integrated over velocity. 

Replacing the evaluated integral in the normalization 
condition gives 



J" oils 



r 2 dr 



o \fi^j 



1. 



(B12) 



The remaining integral can be evaluated quite easily to give 
the definition of : 



16Jbfi s rf /2 



»v 1: - • (B13) 



15^ ' 16r s 5/2 n 

Thus, the complete expression for the DF is 



f(r,v r ) = 



Wrs /2 Q. 



-8 [r s — r — k (v r — Vs) ] 



(B14) 



where we have suppressed the trivial delta functions on vq 
and Vj,. 
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APPENDIX C: DERIVATION OF THE 
PHASE-SPACE DENSITY; VERIFICATION OF 
THE DENSITY DISTRIBUTION 

The phase-space density is obtained by averaging over an 
ensemble of cold distributions, each with a slightly different 
caustic radius: 



f{r,v r ) oc 



exp 



25? 



S \r s + x — r — K (v r — v s ) 2 ] x = X. 

(CI) 



To perform the integral we change variables to the quantity 
inside the delta function, 



u = r s + x — r — K,(v r — v s ) 



(C2) 



which has the same endpoints and a Jacobian of 1 (du = dx). 
Then the integral becomes trivial: 



X = 



so that 



21 



exp 



X = exp 



[it — r s + r + n(v r — v s ) 2 ] 
25? 



[r s — r — k(v t — Vs) 2 ]' 
25? 



5(u)du, 
(C3) 

(C4) 



We can verify that this expression for the phase den- 
sity gives the same functional form as Equation (1) for the 
density p by integrating over v r , although the normalization 
will be different since the phase density is normalized to f . 
The density will be proportional to the integral of I over v r - 

[r s — r — n(v r — v s ) 2 ] 2 



p oc 



exp 



dv r 



(C5) 



We change variables to make the exponent Gaussian in the 
variable of integration, 

_ r s — r — n{v r — v s ) 2 



V2S r 



and determine the Jacobian: 



dv r 



2«(r s — r — y/25 r 



(C6) 



(C7) 



where the upper sign is for the part of the integral where 
v r > v s (the upper half of the parabola) and the lower sign is 
for the part when v r < v 3 (the lower half). This implies that 



we must break the integral into two parts with the following 
endpoints in u, taking into account the sign of the Jacobian: 



-OO < V r < V s 
V a ^ V T < OO 



-oo < u < 



\/25 r 



\/28 r 



So the two parts of the integral are identical; thus we can 

say 



p oc 



(r s -r)/V2S r 



e~ u du 
r s — r — \f28 r u 



(C8) 



This integral can be evaluated as a combination of Bessel 
functions depending on the sign of r a — r. Behind the caustic, 
this quantity is positive; in front it is negative, so we get a 
piecewise solution with the same domains as Equation (1). 
Performing the integral gives the result: 



p oc 



\/\r s - r\e 



4.6'i 



1 B 



(r 3 - r f 



4<5, 2 



(C9) 



which is the same functional form as Equation (1). 



APPENDIX D: EQUATIONS FOR A ROTATED 
CONE 

A cone of height h and opening angle a, with its point at 
the origin and its axis of symmetry along the z axis has the 
parametric equations (in Cartesian coordinates) 



x — (h — it) cos tan a 
y — (h — u) sin 1? tan a 
z — h — u 



(Dl) 
(D2) 
(D3) 



where the parameter u £ [0, h] describes the distance from 
the base of the cone and the parameter £ [0, 27r] describes 
the azimuthal location on the cone. If the cone is rotated so 
that its axis of symmetry points along the unit vector 



h = sin 8 S 



9s cos 4>sX + sin^s sin <j) s y + cos 8 s z 



(D4) 
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then the equations become 

x — {h — u) {cosi? [l — cos 2 4>s (1 — cosfls)] tana 

— sin i? [cos 4> s sin </> s tan a ( 1 — cos 8 S )] 

+ sinfls cos(f> s } 
y — (h — u) {— cos i? [cos (f> s sin S tan a (1 — cos 8 3 )] 

+ sin ■& [l — sin 2 </> s (1 — cos S )] tan a 

+ sin # 3 cos 4> s } 
z — (h — u) [— cos ($ — i/) s ) sin# s tana + cos# s ] 



where I? is the discriminant 



(D5) 



(D6) 
(D7) 



To determine z m i n and z m ax as a function of x and y, we solve 
the system of the x and y parametric equations to obtain 
the parameters u and The system can have zero, one, 
or two solutions depending upon the values of x and y: zero 
solutions for lines of sight that do not intersect the cone, one 
solution for lines of sight that intersect the cone once, and 
two solutions for lines of sight that intersect the cone twice. 
In cases where there is one solution, the other limit can safely 
be taken to be ±oo, where the plus sign is taken if the bulk of 
the cone is in front of the intersection point (the solution is 
the lower bound of z), and the minus sign is taken if the cone 
is behind the intersection (the solution is the upper bound 
of z). The system can always be solved analytically and we 
present the solution here. First we will define some auxiliary 
quantities to make the notation simpler. We extract the 9 S , 
<f> 3 , and a dependence of Equations (D5-D6) into coefficients 
that need only be calculated once for a given cone: 



A x 
B x 

Cx 

Ay 



tan a [l — cos 2 <f> s (1 — cos 6 S )] 
tan a cos (j> s sin 4> s (1 — cos 8 S ) 
sin 6 S cos d> s 



— B x 



= tan a [l — sin 2 <f> s (1 — cos 9 S )] 



Cy 



sm a s sin t 



so that Equations (D5-D6) become 

x — (h — u) (A x cosi9 — B x sini9 + C x ) 



y 



(h - u) (-A y cos i9 + By sin ■& + C y ) . 



(D8) 
(D9) 
(D10) 
(Dll) 
(D12) 
(D13) 

(D14) 
(D15) 



We can solve for ■d by dividing the two equations, noting 
that if x = or y = then we are at the point of the cone 
(u = h) and § is degenerate. We obtain an equation for ■& in 
terms of the ratio r\ = x/y: 



5sini9-Ccosi9 + /C = 0, 



where 



C — A x 4~ ^}A y 

S = B X + TjBy 

K, = — C X +7]Cy. 



(D16) 

(D17) 
(D18) 
(D19) 



Equation (D16) can be recast as a quadratic equation 
for either sin$ or cosi9. To avoid using inverse trigonomet- 
ric functions (and the associated difficulties in choosing the 
right branch) we simply solve for both and use them in the 
rest of the solution: 

(cos*) ± = ^±4? (D20) 



(sintf) ± = 



C 2 +S 2 
K,S±CV 
C 2 +S 2 



(D21) 



V 2 = C 2 +S 2 -K 2 . 



(D22) 



As usual, if T> < the point is outside the cone and the 
equation has no real roots; otherwise it has two real roots. 
When two solutions exist for •d, one or both of them may 
lead to a value of u outside its allowed range. 

Having determined # either Equation (D5) or Equation 
(D6) can be used to determine the quantity h — u that is 
necessary to find the limits on z: 



[h-u) i 



A x (costf) ± - B x (sintf) ± + C x 

V 

-Ay (costf), + By (sin??), + C y 



(D23) 



The quantity h — u should be in the range (0, h] — if it is not, 
that value of $ is discarded as a root and one of the limits 
in z goes to ±oo. 

Finally, the limits z± are given by plugging in the valid 
solutions for u and 

z± — (h — u) ± {cos6< s — sin^s tana [(cosi?) ± cos<^> s 

+ (sintf) ± sin<^]} (D24) 

By inspection, we see that z m i n is not always equal to z- and 
2max is not always z + since both (cosi?) ± and (sin$) ± can 
take any sign: the roots must be compared and the smaller 
assigned to z m i n . If one root is out of range in u, its out- 
of-range z value can still be calculated (for this work the 
height h of the cone is arbitrary because the density function 
effectively cuts off the cone along a spherical segment) and 
compared to the value of the in-range root to determine 
whether the out-of-range limit is z m i n or z max . 

In constructing a fitting routine, the partial derivatives 
of E with respect to the projection parameters a, S , <f> a and 
to the density profile parameters r s ,8 r , k, fo may be needed. 
In the following, we define the shorthand 



[C0S(?? — (f> s 

[sin(# — (f> s 



( COS la, COS fa + ( Sin la, Sm 4>s 
( Sm ^)max COS <t>« ~ ( COS ^)ma X Sm ^ 

are taken to be the roots 



where (cosi9) max and (sintf) max 
that lead to the value of z max , and likewise 



[cos($ — 4> s 
[sin(i? — (f> s 



(cos i9) min cos 4> s + (sin •i?) min sin 4> s 

( Sin ^)min COS 08 - ( COS ^)min Sm ^» 



for the roots leading to the value of z m i n . A similar notation 
is used to denote the appropriate root of (h — u). We also 
use the shorthand 



p(z max ) = p(\Jx 2 + y 2 + z2 iax ) 



(D25) 



in the following, since x and y are understood to be constant. 
With these definitions, the derivatives with respect to the 
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projection parameters are 
<9E 



{p(Zmin) [C0S(;& - 0s)] min (ft - u) n 



da T 

-p(Zmax) [cOs($ - 0s)] max (ft - u)max} (D26) 
= Y {p(2min) [(ft - w)min Sill S 

+ [cos($ — s )] min cos# s tana] 
-p(zmsx) [(ft - u) max sin 6» s 

+ [cos(tf - s )] max cos 6» s tana]} (D27) 
<9E sin(9 s tana f . 

= Y |P(^min)(ft - Wjmin [Sin(l? - s )j min 

-p(Zmax)(fe - «)mix [sin(l? - 0s)] max } (D28) 

The derivatives with respect to the profile parameters 
are all integrals of the derivatives of the density, of the form 

o7Vi T J Zmia 8-Ki 

for a given parameter 7T; , since the limits of the integral and 
the variable being integrated over do not depend on any of 
the 7r; and the function p, although defined as a piecewise 
function in Equation (1), is continuous over the entire in- 
tegration range. To compactly write derivatives of p with 
respect to the parameters, we expand the definition of B to 
include other Bessel functions: 

R r „N = / f [Z(2n+l)/4(«) +X_(2„+l)/4(u)] T < T s 

tt nW - | ( _ ir+ i f [ X(2ri+1)/4 ( U ) -I_ (2 „ +1)/4 ( W) ] r > rs 

(D30) 

with u = (r — r s ) 2 /4(5^ as before. The definition used in 
Equation (1) is equivalent to Bo in this new notation. With 
this simplification, the derivatives dp/dm are: 

dp _ fo c -u V\r-r s \ x 



dr a V27TK 

x |u [2Bb(«) - Ba(u) - Bi(u)] - ^Bol ») }■ (D:iJ ) 



<9<5 r V27TK '' <5r 

x [2B (it) - £ 2 (u) - Si (it)] (D32) 

1 - -k « D «» 
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